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In the study of complex networks (systems), the scaling phenomenon of flow fluctuations refers 
to a certain power-law between the mean flux (activity) (Fi) of the ith node and its variance 
<Ji as Gi oc (Fi) a . Such scaling laws are found to be prevalent both in natural and man-made 
network systems, but our understanding of their origins still remains limited. In this paper, a non- 
stationary Poisson process model is proposed to give an analytical explanation of the non-universal 
scaling phenomenon: the exponent a varies between 1/2 and 1 depending on the size of sampling 
time window and the relative strength of the external/internal driven forces of the systems. The 
crossover behavior and the relation of fluctuation scaling with pseudo long range dependence are 
also accounted for by the model. Numerical experiments show that the proposed model can recover 
the multi-scaling phenomenon. 



PACS numbers: 



I. INTRODUCTION 



The studies of the complex networks have attracted in- 
creasing interests since the last decade tiJ» [2)> [2)> ]£. Much 
work has been devoted to the understanding of the dy- 
namic processes of network flows and significant advances 
has been achieved. However, there are still many ques- 
tions which remain to be answered thoroughly. Among 
them the origin of the fluctuations of network flow is of 
particular interest, which will be addressed in this paper. 

The study of the scaling of fluctuations can be dated 
back to the Taylor's law initialized in [5j]. Recently, 
related studies attract increasing interests again, since 
Menezes and Barabasi's proposal of the so called scaling 
analysis. The recent research aims to describe the dy- 
namics of a large number of nodes simultaneously and ex- 
amine the collective behaviors of the systems ^' [3> 2 . 
They investigated the coupling between the average flux 
and fluctuations and found that in complex networks 
there exists a characteristic coupling between the mean 
and variance of flux on individual nodes as oi oc (Fi) a . 
Such scaling law was found to hold in diverse natural 
and man-made systems. Moreover, they claimed that 
real systems belong to one of two distinct universality 
classes. The first class yields a ~ 1/2, which indicates 
that system dynamics are dominated by internal factors. 
Internet and microchip systems belong to this class. The 
other class yields a w 1, which indicates that the fluctua- 
tions mainly come from external driven forces. Highway 
traffic, river networks, and the World Wide Web belong 
to this class. 

These interesting findings are soon followed by fur- 
ther studies. The scaling phenomenon has been con- 
firmed in a wide range of complex networks, e.g. stock 
markets 112). liiJ, gene systems lis, download networks 
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li2J, software developments and urban transporta- 
tion networks 112). In these systems, a rich variety of 
exponent a from 1/2 to 1 can be observed, which ap- 
proach 1 as the size of sampling time window or the rela- 
tive strength of the external driven forces to the internal 
dynamics increases. Some systems even exhibit the so 
called multi-scaling behavior (i.e., the power-law extends 
to the qth moments of Fi) and the crossover (i.e., dif- 
ferent groups/types of nodes in the same system possess 
different values of a) behavior [i2J. lilJ. li2J. liSJ. These 
phenomena are believed to reflect the competition be- 
tween the external and internal dynamics, as well as the 
inhomogeneity of the systems' structure [13 ■ LIS ' [13. 

Several models have been proposed to explain the 
underlying mechanisms that generate the scaling laws 
[3,[8],[i6j,[ir],[i8],[i9] i most f w hich are based on ran- 
dom walk or diffusive dynamics models. These mod- 
els can reproduce the scaling phenomenon and link the 
dynamic processes of networks with the associated net- 
work topologies. Usually, the following questions are ad- 
dressed: 

1 . The scaling law is better to be represented via both 
explicit analytical expression and numerical simu- 
lations E). [13. [IS]. [IS. [20]. [21]. 

2. The non-universality values of a from 1/2 to 1, the 
effect of the sampling time windows and the inter- 
nal/external driven forces should be explained in 
an integrated framework 2)' 113 ' [23 . 

3. The crossover phenomenon should be accounted for 

[9], [20] [21] _ 

However, we are still interested in finding a even more 
simpler model other than random walk/diffusive dynam- 
ics models to answer the above questions. In the authors' 
earlier attempt, a non-stationary Poisson process model 
was proposed in [2(| to overcome the above shortcom- 
ings. This model could explain, both analytically and 
numerically, the transition of a from 1/2 to 1, as well as 
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its dependence on the sampling time windows and exter- 
nal/internal driven forces. A similar model also appeared 
in [2l[. However, both [2(J and [2l[ assumed that the 
length of the intervals between two consecutive jumps of 
the arriving rate equals that of the sampling time win- 
dow. In order to better model the phenomena in real 
systems, it is necessary to relax this restriction and allow 
different size of the jumping intervals and the sampling 
time windows. 

This paper extends the model in [2(| to more general 
cases and provides a more comprehensive study. In par- 
ticular, analytical solutions are provided for a's depen- 
dence on sampling time windows and external/internal 
driven forces, as well as for the crossover phenomenon. 
Moreover, the connection of fluctuation scaling with 
pseudo long range dependence is explained. The multi- 
scaling behaviors are also investigated numerically. The 
proposed model is able to capture the essential mech- 
anism that drives the fluctuations of flow and offer a 
unified and concise picture for the various phenomena 
observed in real traffic flows. 



II. THE NON-STATIONARY POISSON 
PROCESS MODEL 

Suppose the network system consists of N nodes, which 
will be indexed as i = 1, N in the rest of this paper. 
The arriving flow fi(t) of the ith node at the tth time 
interval is assumed to behave as a non-stationary Poisson 
process 

Pr(Mt)=n) = e - (1) 

where n = 1, 2, ... . Pr denotes the abbreviation of the 
probability. 

Moreover, the varying process of the arriving rate Aj(i) 
is assumed to change according to a network-wide finite- 
state deterministic or stochastic process. The system will 
enter one state during a certain time interval r and then 
enter another state, which as a result changes the arriving 
rate as well as the arriving flux of the Poisson process. 

For instance, the widely used Markov-modulated Pois- 
son process (MMPP) model is a typical example of such 
models. It is a doubly stochastic process where the in- 
tensity of a Poisson process is defined by the state of 
a Markov chain. The Markov chain can therefore be 
said to modulate the Poisson process, and thus comes the 
name MMPP in many literatures [22). E2J. [2£] [25] [26], [27] _ 

Fig. I shows a Poisson process modulated by a two-state 
Markov process, in which we can observe the shift of av- 
erage arriving flow every r. 
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Fig. 1. Examples of a Poisson process modulated by a 
two-state Markov process. 

Before discussing the characteristics of such models, 
let's introduce some phrases to distinguish the different 
time intervals employed in the rest of this paper. 

1. Intrinsic time interval denotes the minimum dis- 
cernable (from the viewpoint of observer) time in- 
terval in the model, which is determined only by 
the non-stationary Poisson process. And we use 
the integer number set {t = 1,2,3, ...,T} to index 
the intrinsic time points. 

2. Endogenous jump time interval denotes the last- 
ing time after the process enters one state and be- 
fore it leaves it, which is determined only by the 
non-stationary Poisson process, too. Here, we as- 
sume it to be an integer multiple of the intrin- 
sic time interval and use the integer number set 
{j = 1, 2, 3, T/t} to index the jump time points. 

3. Exogenous sampling time interval denotes the 
length of the sampling window (from the view- 
point of observer), which can be changed by the 
observer. We assume it to be an integer multiple 
of the intrinsic time interval, too. We use the sym- 
bol At to represent it and the integer number set 
{s = 1, 2, 3, T/At} to index the sampling time 
points. 

Suppose the dynamic varying process of the arriving 
rate Xi(t) can be depicted as 

Ai(t) = kiX(j), for (j-l)r + l<t<jT (2) 

where ki is a scale coefficient which accounts for the fac- 
tors controlling the relative magnitude of flow on node i. 
For example, in a random walkers model, fcj stands for 
the degree of node i [2iJ . 

A(j) represents the network-wide stochastic scalar pa- 
rameter at the jth jumping interval. We assume all the 
nodes in the network systems follow a synchronized tran- 
sition tempo so as to reproduce the network- wide fluctua- 
tions in a simple way. Studies show that we can relax this 
assumption to reproduce even more complex phenomena. 
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In this paper, we mainly focus on the steadily state of 
this process, since it determines the steady distribution of 
the scalar parameter X(j) and thus shapes the mean and 
variance of flux fi(t) at the ith node. Suppose we have 
altogether M possible values of the scalar parameter X(j), 
which are denoted as {\ m \m — 1,...,M}; and each A m 
corresponds to a certain state of this process. Moreover, 
we assume A(j) have the following stationary distribution 



Pr (A(j) = \n) = Pn 



(3) 



In practice, we usually deal with sampled flows series 
instead of the original series /, (t) . Given a sampling win- 
dow of length At, the observed flow data at the sth sam- 
pling index is the summation of the original flow values 
within the sampling window: 



sAt 

E 

t=(s-l)At+l 



fi(t), for s = l,...,T/At (4) 



Under the assumption of crgodicity, the time average 
(■} of the flow series equals its equilibrium value. Thus 
the observed average flow at node i can be written as 



T/At 



T/At sAt 



r/ ^E^ At (*) = ^E. X. m 



= 1 t=(s-l)At+l 



= t&(*) = tE E /*(*) 



t=l 

T/r 



T 

At 
~T 

3 = 1 

hAt(X) 



3=1 t=(j-i)r+i 

T/r 



J2k i X(j)r = k i At^j-J2x(j) 

3 = 1 



(5) 



where (A) := ^ Y,Jlt Hi) = Em=i x ™p m - 

To derive the variance (c^*) 2 of the observed flow, we 
need to distinguish two cases in terms of sampling win- 
dow size At. 



At 



Fig. 2. Diagram of Case I), where At < r. 



Time 



I) At < t (see Fig. 2). Further assume that t = nAt, 
k £ N, thus the arriving rate keeps constant in each sam- 
pling window. Noting that the sth observed data satisfies 



F t At (s) ~ P(kiXjAt) for (J - 1)« + 1 < s < jn, we have 



(^At\2 

l°i ) 



T/At 

^E [*i At (*)-<*i At >f 



At 



T/r j K 



^E E -<*s 



At\2 



3 = 1 s=0'-l)*+l 



T 
Ai 



~ E ^ 



{hAt(X)Y 



3 = 1 
T/r 

K [hX(j)At + (k t X(j)At) 2 ] - (hAt(X)) 2 



s=(j-l)«+l 



T 



3 = 1 



= k t At{X) + kf(At) 2 [VarA + (A) 2 ] - (fc 4 Ai(A)f 

= kiAt(X) +fc 2 (At) 2 VarA 
VarA 

W 



(F^ + ^iF^r 



(6) 



where VarA := f Ejii M) - (A)]* - E™=i A^p 



M x2 



2^m=l A mPn 



At 



Time 



Fig. 3. Diagram of Case II), where At > r. 

II) At > t (see Fig. 3). Further assume that Ai = kt, 
k 6 N \ {1}. We have 

T/At 



At 



At\2 



s=l 
T/At 



At x ^ 



s=l 



E E /*(*) 

j=(a-l)K+l i=(j-l) T +l 



(7) 

Let A s := i X^=(s-i) K +i Mi) denote the average of 
A(j) within the sth sampling window. Noting that 
(A s ) = A* £j/f * A s = (A), Eq.0 can be simplified as 

= fciAt(A s ) + fcf(Ai) 2 [ K VarA s + <A S } 2 ] - {k t At{K s )f 

(8) 



Moreover, we have 



T/At 

VarA s = r/At E ( A * ~ < A *» 2 = ~ VarA ( 9 ) 

' 6 — 1 
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Thus, we can summarize Eq.©, ([6]) and (JH) as 

(F^) = kiAt(X) (10) 

and 

{at r = hAt(X) + (fc 4 Ai) 2 VarA = + |^ At } 2 

(11) 

Eq. (jllj) indicates that no matter the comparative ratio 
between the size of sampling window and the size of the 
jump interval, the variance (c^*) 2 is always a compound 
of the mean F^ 1 and its power of 2. Therefore, we can 
always find the scaling exponent a varies between 1/2 
and 1. 

Besides, it should be pointed out that although we 
assume k to take integer values, a non-integer k does not 
affect the main conclusions of the proposed model. 

III. THE SCALING OF FLUCTUATIONS AND 

PSEUDO LONG RANGE DEPENDENCE 
REPRODUCED BY THE PROPOSED MODEL 

Based on Eq. (TIT)]) and (jTTJ) , the scaling law of the pro- 
posed model can be analytically explained. Particularly, 
we will discuss the three scaling phenomena which have 
been frequently observed both in simulations and real 
systems in previous literature. The connection to pseudo 
long range dependence will also be investigated. 

1) The effect of the external and internal forces 

The time-variation of X(j) can be regarded as the re- 
sult of external driven force; thus the ratio VarA/(A) re- 
flects its relative strength to internal dynamic factors. If 
the external driven force is weak, we have VarA/(A) <C 
l/(fciAi); thus the second term in the r.h.s. of (fTTj) can be 
omitted, yielding a = 1/2. On the other hand, if the ex- 
ternal driven force is strong and VarA/(A) 3> l/(kiAt), 
the first term can be omitted, yielding a = 1. For in- 
termediate values of VarA/(A), the scaling law actually 
breaks down 122], but we still have effective values of a 
ranging between 1/2 and 1. 

Particularly, we are interested in network sys- 
tems containing two or more driven forces (some are 
global/system- wide and external; some are local/district- 
wide) . To separate these might help distinguishing which 
effect is dominant at a certain time, which will further as- 
sist us to take appropriate actions/controls. For example, 
in a previous report tiSJ , we found the scaling phenom- 
ena in urban ground traffic network might be explained 
as the interaction between the nodes' internal dynamics 
(i.e. queuing at intersections, car-following in driving) 
and the changes of the external (network-wide) traffic 
demand (i.e. the every day increase of traffic amount 
during peak hours and shocking caused by traffic acci- 
dents). This may allow us to further understand the 
mechanisms governing the transportation system's col- 
lective behavior. 

2) The influence of the sampling window 



If the size of the sampling window At <C (A)/(fciVarA), 
we can omit the second term in the right-hand side of 
(|11[) . which contain the quadratic form of At, and thus 
a = 1/2 [34|. If At > (A)/(fc l VarA), we can omit the 
first term in the right-hand side of (jll[) and thus a = 1. 
For intermediate values of At, the effective a lies between 
1/2 and 1. 

3) The crossover behavior 

For nodes with small ki, the second term in the right- 
hand side of (Till can be omitted and it yields a = 1/2. 
For nodes with large ki, the first term in the right-hand 
side of (jTTJ) can be omitted and thus a = 1. If the system 
has heterogeneous structure and consists of a broad range 
of ki across the nodes, it can be expected that different 
subsystems (i.e., different groups of nodes with different 
ki and thus different (F^)) have different values of ex- 
ponent a. 

4) The coexistence with pseudo long range dependence 
Besides the scaling law of fluctuations, the long range 

dependence represents another kind of power-law, which 
couples the variance and the size of the sampling window. 

a t {At) cx At H{ ^ (12) 

where H(i) denotes the well-known Hurst exponent. 

The coexistence of these two power-laws has received 
considerable interest, since they describe the behavior 
of the same standard deviation <j,(At). Previous stud- 
ies showed that these two power-laws coexisted in some 
systems (e.g., stock markets), and several models were 
developed to explain such phenomenon tiZJ> 122]' 122]. 

In the proposed model, the power-law in (TT2"]) can be 
derived directly from (|12p . with H(i) ranging between 
1/2 and 1, which recovers the possible coexistence of 
long range dependence (actually "pseudo long range de- 
pendence" here) and the scaling of fluctuations. More- 
over, ([TT]) predicts that H(i) -v 1/2 when VarA/(A) -> 0, 
which reflects the uncorrelated nature of a simple Pois- 
son process. On the other hand, when VarA/(A) > 0, the 
time- variation of the arriving rate induces long-range de- 
pendency to the modulated Poisson process and results 
in H{i) > 1/2. 

It should be pointed out that the long-range depen- 
dency discussed above is indeed the so called "pseudo" 
long range dependency (PLRD) instead of mathemati- 
cally rigorous LRD. Usually, LRD means that the decay 
of the autocorrelation function is hyperbolic and decays 
slower than exponentially for all the time lag (time scale) . 
However, long range dependency properties, heavy tail 
distributions, and all other characteristics of real time 
series (i.e. Internet traffic, ground vehicular traffic flow) 
are meaningful only over a limited range of time scales. 
Many recent approaches 122J 123 12£] had illustrated that 
appropriately constructed Markov models appear to be a 
viable modeling tool in the context of modeling LRD time 
series over several time scales. Particularly, the MMPP 
models are shown to provide good matches of LRD prop- 
erties under large time scales. The key idea behind is ob- 
taining a model that has fast transients on a local scale 
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and evolves with a slower time constant between differ- 
ent groups of states. In other words a model that has 
some implicit form of memory, which is the origins of 
such PLRD. The model proposed here yields some kind 
of PLRD also via this way. 

Besides, the proposed model also indicates clear non- 
universality, which means that the Hurst exponents of 
different nodes may vary even though the underlying 
Poisson mechanisms for these nodes are the same. Thus 
we should be careful when applying concepts like scaling 
and universality to complex systems 122) . 

IV. NUMERICAL EXPERIMENTS 

In this section, some numerical experiments are de- 
signed to verify the proposed model. 
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Fig. 4. The transition behavior of a w.r.t the compe- 
tition of the external/internal driven forces produced by 
simulation. 

In order to observe the influence of the strength of 
the external and internal forces, the corresponding pa- 
rameters of the model are set as T = 100000, r = 1, 
At = 1. N = 20, k = [1,2,3, ...,20]. M = 10, the ra- 
tio VarXj (A) is controlled to vary from 0.01 to 200. In 
the first test. The transition of a from 1/2 to 1 under 
different VarA/(A) can be clearly observed from Fig. 4. 
Similar results can be found that At > r cases and are 
thus omitted here. 
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Fig. 5. The transition behavior of a w.r.t the sampling 
time window produced by simulation. 



In the second test, the corresponding parameters of 
the model are set as T = 1000000, r = 1. N = 20, 
k = [10, ...,20]. M = 2, p m = 1/2 for m = 1, 2, A m = 
[0.8, 1.1]. The sampling window size At is controlled to 
vary from 1 to 15. As shown in Fig. 5, a increases as At 
is increased. 




<Ff< > 



Fig. 6. The crossover behavior produced by simulation. 



To test the existence of crossover behavior, the cor- 
responding parameters of the model are set as T = 
100000, r = 1, At = 1. M = 10, p m = 1/10 for 
m = 1, 10, A m = [1,2, 3, ...,10]. N = 20, k = 
[0.5, 0.5 2 , 0.5 8 , 1, 2, 5, 5 2 , 5 10 ] in the third test. The 
result is plotted in Fig. 6, which shows clear crossover 
behavior. 

Another important phenomenon numerically repro- 
duced by our model is the multi-scaling behavior, which 
is the extension of the original scaling law to the gth- 
order central moments 



^ = (l^-(^>IV(0 ,a(,) (13) 



Generally, this multi-scaling law (i.e. a dependence 
of a(q) on q) is assumed to be related with the multi- 
fractality of time series liZJ. In the fourth test, the corre- 
sponding parameters of the model are set as T = 10000, 
t = 10. M = 10, pj = 1/10 for j = 1, 10, 
Xj = [2, 3, ...,11]. N = 20, k = [1,2,. ..,10]. Fig. 7 plots 
the three q-a(q) plots for At = 1, At = 10 and At = 100, 
respectively, in which a{q) heavily depends on q. 
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1.4 




ations in complex systems. The scaling exponent a with 
non-universality values between 1/2 and 1 can be analyt- 
ically derived. The influences of the sampling time win- 
dow and the external/internal force ratio, the crossover 
behavior, the multi-scaling phenomenon, and the con- 
nection with long range dependency are also naturally 
explained. The model is verified by numerical simula- 
tions. 



-20 -15 -10 -5 5 10 15 20 

q 

Fig. 7. The multi-scaling behavior produced by simula- 
tion. 



V. CONCLUSION 

In summary, a non-stationary Poisson process model is 
introduced in this paper to explain the scaling of fluctu- 



Our work also sheds light on the well-known debate on 
whether Poisson process models are applicable in model- 
ing complex network systems which often exhibit bursti- 
ness and long range dependency sU' 122J' [22J. The results 
presented in this paper support the argument that the 
non-stationary Poisson process are still powerful in mod- 
eling the multi-scale behaviors of complex time-variant 
systems with complicated interactions between external 
and internal dynamics. 
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